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A Method for Finding Roots of Arbitrary Matrices 


1. Jacobi’s Method for Real Symmetric Matrices. There is a well known 
method due to Jacosr' for diagonalizing real symmetric matrices. It consists of 
performing a sequence of orthogonal transformations (rotations), each one on a 
two-dimensional subspace of the underlying vector space in which the matrix is 
considered to be defined. Thus, if we call the original matrix A, we transform with 
an orthogonal matrix Sam) as follows: 


(1) B= STemyA Stem): 


(The superscript T indicates the transposed matrix.) The main idea of Jacobi’s 
method is to annihilate one of the off-diagonal elements (A:,) by this rotation. 
After we have annihilated A;,, with this rotation, we go to element Ax», and find 
a transformation Sm) which will annihilate A,-,,. This last transformation, if 
it affects Aim, will in general ‘‘deannihilate”’ it. Nevertheless, one selects each off- 
diagonal element in turn and performs a transformation which will annihilate it. 
After one has covered every off-diagonal element once, one has performed what 
we will call an iteration. The essence of Jacobi’s method is that even though one 
undoes the annihilation (in general) of all the off-diagonal elements except the 
last one, with a sufficient number of iterations the off-diagonal elements will 
§ converge to zero, leaving a diagonal matrix. One obtains the final transformation 
matrix by post-multiplying by the San) at each annihilation, starting with the 
unit matrix. Thus: 


(2) Stinat = *** SkemeStom *** Skeme 


where the superscript refers to the number of the iteration and R is equal to the 
number of off-diagonal elements (= $n(m — 1) for an mth order matrix). This S 
has the property that: 


(3) STAS = D 


where D is the diagonal matrix whose elements are the roots of the matrix. With 
each root, there corresponds an eigenvector, which is the column of the S-matrix 
which occupies the same position in the S-matrix as the column containing the 
eigenvalue in question occupies in the D matrix. 

2. Generalization for Arbitrary Matrices. There are various fundamental 
theorems which show that a real symmetric (or, for that matter, Hermitian) 
matrix has real roots and can always be diagonalized by an orthogonal trans- 
formation of some kind. 

These theorems all break down for non-Hermitian matrices, real or otherwise. 
In fact, it is known that there exist matrices which no collineatory transformation 
can diagonalize. These are all degenerate, and the best reduction that can be 
effected is to the Jordan canonical form. However, if we give up the ideal of 
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diagonalizing matrices, and restrict ourselves to triangularizing them, then we 
have at our disposal a theorem of ScHurR? which states that an arbitrary matrix 
can be triangularized by unitary transformations. Therefore, instead of specifying 
Stam) to be orthogonal, we shall require it to be merely unitary, and perform the 
following transformation on the arbitrary complex matrix, A: 


(4) B = 


The elements of Sum) are now complex, and it is in general impossible to 
characterize S by one parameter. There are really two parameters necessary to 
characterize the transformation. We shall define the 2 X 2 non-trivial submatrix 
of Sam) as follows: 


Sek Sim | | @ —€ 
where a is real and positive, and the following relation holds: 


(6) a? + |c|? = 1. 


It can be shown that there is no loss in generality in thus restricting Sam). 
The elements of the transformed matrix are given by: 


Bu = aAu + CAim 
Bim = — + 


(7) 


(8) Biz = + CAnj 

CAnj GA mj 

Bee = + |C|*Amm + ACA bm + mx 
(9) Bim = — + — Arr) 


Buz = — + 6(Amm Arr) 
Bum = @Amm + |c|*Ae — acAnm — 


In view of the non-hermiticity of A it is in general not possible to make both 
Bmx and Bim vanish. However, it is always possible to make B,, alone vanish, 
since there are three equations in the three unknowns a, Ri(c), Im(c) if we bear 
in mind that setting B,,. = 0 gives two equations and that (6) must be satisfied. 


By setting: 
(10a) a=[1+ 
(10b) c= ali + = 


we can reduce the number of unknowns and equations to two. We then have, for pu: 
(11) Amp? — (Amm — Arr) — Am = 0; 
hence (setting Amm — Ar = 2A4nx), we obtain: 


(12) = VAR, + me: 
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Of these two roots, we choose the one of smaller modulus, in order to avoid large 
rotations when possible. (This can be done without much trouble in a digital 
program.) 

There are properties of these unitary transformations which are of impor- 
tance, viz.: 


(13) |Bual? + |Bim|? = + |Aim|? 
(14) + | = + 
(15) | Brel? + |Bem|? + |Bme|? + | Bm |? 
| + + | Ame]? + | 


It is no longer possible, as in the real symmetric case, to prove that the sum 
of squares of absolute values (S.S.A.V.) of the off-diagonal elements decreases. 
What is necessary to prove is that the S.S.A.V. of the sub-diagonal elements 
decreases to zero. It is not even true that this sum decreases monotonically.? 

An example to illustrate this is the following: 


(16) Az=]01 1|. 
One may annihilate the lower left element with the transformation : 
(17) Say = 010 


0 0 
but the result is: 


(18) SévA Se =| 


and we see that the S.S.A.V. of the sub-diagonal elements has increased from 
1 to 2. However, one simply keeps on transforming away sub-diagonal elements 
without regard to the fluctuations of the S.S.A.V. of these elements. Finally they 
should decrease to zero and the process converge. 

The justification for the last statement is based on “experimental” evidence, 
gained by triangularizing a variety of matrices with a program written in the 
“Speedcoding”’ floating-point system for the IBM 701. Floating-point is unneces- 
sary in this problem and slows the machine by a factor of more than ten, but this 
program was written originally as much to test Speedcoding as to test triangulari- 
zation. The average time per multiplication or addition in this system is 4 
milliseconds. 

With a series of matrices’ whose elements (both real and imaginary parts) 
were picked at random, the following are the number of iterations and times re- 
quired for convergence to 10 decimal places (actually, with floating point, to 
10-*° for sub-diagonal elements). 

These data, though incomplete, indicate that the number of iterations in all 
probability does not increase more strongly than the first power of the order of 
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TABLE I. Time for Convergence 


Order Time per Iteration Number of Iterations 
2 2 3 
3 5 6 
10 9 
5 20 12 
6 40 14 
7 55 14 
8 85 17 


the matrix. For degenerate matrices, the number of iterations increases consider- 
ably, rising to about 25 for a fourth order matrix with four equal roots. It is 
noteworthy, also, that the roots (which are, of course, the diagonal elements when 
triangularization is complete) are not accurate. This is because the unitary 
transformation really preserves the invariants of the matrix, i.e., the coefficients 
in its characteristic equation, and these do not specify the roots uniquely because 
of rounding error. For example, let 4) be a double root of : 


(19) (A — = — + Ac? = O. 


What we are given originally, in general, is the expanded form of this equation. 
Now let us substitute for (Ao, Ao) the pair (Ao + €, Ao — €). The characteristic 
equation is then: 


(20) (A — Ao — e)(A — Ao +) = — + — =O. 


Now if |\o| ~ 1 and |e| ~ 10-5, then the characteristic equation in (20) will 
differ from that in (19) by é&, or a quantity ~10-"° which is of the same order of 
magnitude as the rounding error in a ten-digit machine. Double precision opera- 
tions will not increase the accuracy, since the coefficients of the characteristic 
equation are not known with more than single accuracy. The roots, as we see, 
may be specified with an ambiguity much greater than the least count of the 
machine. In general, if there are m equal roots and D digits, the roots will be 
accurate to about D/m digits. Therefore, it is in general impossible to decide 
numerically whether roots are equal or just very closely spaced. 

These remarks do not apply to Hermitian matrices, where equal roots will 
come out equal within the actual rounding error of the machine. There is evidently 
some connection between the ‘‘degree of Hermiticity’’ of a matrix and the accu- 
racy with which its repeated roots can be found. 

It is of interest at this point to indicate the complications that arise in at- 
tempting to find the eigenvectors of the triangularized matrix, which distinguish 
markedly the Hermitian and non-Hermitian cases. In particular, it will be shown 
that, when a root is repeated, it is sometimes impossible to obtain the full number 
of eigenvectors corresponding to that root. First we look at the Hermitian case. 

The original eigenvalue equation was: 


(21) = 
which now becomes: 
(22) = a® 
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with 
(23) = 
Since Diy = a5;; we have: 


(24) = ay = 


7 


Hence, ¢: must be equal to 4,;, i.e., all components of ¢ vanish except the 
kth. Hence also, from (23): 


(25) V™ = Sie = Sides = Sa. 
i i 
This shows why y is the kth column of S. 
Now, in the case where A is arbitrary, we can at most triangularize it to obtain: 
(26) St*AS = T 


where TJ is a triangular matrix. Analogously to the Hermitian case, the trans- 
formed eigenvalue equation is: 


which is written out as: 


N 
(28) Tie = a” 


We subtract one side from the other to obtain: 


N 
(29) LX (Ty — = 0. 

In eq. (24), we had in place of T;;, a“8,;;. This gave (a — a) go, = 0. Here, 
however, we have terms other than the diagonal ones coming in, so the situation 
is not so straightforward. Eq. (29) is, in fact, a triangular system of equations for 
¢;” which must be solved for each value of k. Let us take k = N, so that the last 
diagonal element vanishes. Let also i = N and consider the last equation, 
which is: 


(30) (a™ — a) = 0, 


Since the first factor vanishes, gy“ is arbitrary. Hence we may solve the tri- 
angular system for all the other components of g™ in terms of gy“ (assuming 
no degeneracy, i.e., that there is no other vanishing diagonal term). Next we set 
k = N — 1. The last equation of the system (29) is: 


(31) — er» = 0. 


Since the first factor does not vanish (by the assumption of no degeneracy), the 
second factor must. The (N — 1)th equation of (29) is: 
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Here, the second term vanishes, and so does the first factor of the first term. 


Hence ¢{" 5” is arbitrary. We may then solve for the other components of g-” - 
in terms of g{Y >”. Similarly, we may calculate all the vectors ¢ in terms of 8 
arbitrary scale factors, which may be chosen so as to normalize the (transformed) tat 
eigenvectors o™. Finally, from yY™ = Sy, we obtain all the eigenvectors of the thi 
original problem. In 
When we do have degeneracy, let us, for definiteness, place one of the repeated 
roots in the last position of the diagonal. Now we proceed as above, except that 
for some value M (say) of the index k we will come across a diagonal coefficient 
in our triangular system of linear equations which vanishes just as the last diagonal 
element did. Two possibilities are now open. Either the rest of the linear combina- 
tion (exclusive of the diagonal term) which constitutes the Mth equation vanishes, 
or it does not. If the former is true, then we are entitled to choose for gy 
‘another arbitrary number and hence obtain another eigenvector. If the latter is fot 
true, then we are forced to make it equal to zero by setting gy“) = 0. This ie 
causes us to lose one of the basis vectors of the degenerate subspace defined by cel 
the repeated root. We may then start as before by setting gy“ equal to some for 
arbitrary number and solve for the rest of the components of g” as before. The int 
case when we lose one of the eigenvectors belonging to an eigenvalue corresponds m: 
to the case in the Jordan canonical form when a 1 appears attached to two equal 
roots. When there is no 1 attached to a repeated root, we then have the full 
complement of eigenvectors for that root. Similar considerations apply to roots 
of higher multiplicity than 2. 
J. GREENSTADT 
International Business Machines Corporation 
New York, New York 
1See, for example, R. T. GrEGory, MTAC, v. 7, 1953, p. 215. 
2?F,. MuRNAGHAN, Theory of Group Representations. Johns Hopkins Univ. Press. Chap. 1. 
* Unless a suggestion of J. von NEUMANN (private communication) is followed, according to 
which one attempts to minimize the S.S.A.V. of the sub-diagonal elements in the mth row and 
kth column between the pivotal subdiagonal element and the diagonal (including the pivotal 
element itself). This method, however, is much more time-consuming than that proposed in 
this paper. 
ce 
Tl 
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1. Abstract. For making interpolations at different parts of a table, double sh 
interpolation formulae for mixed forward, backward, and central differences have fo 
been derived. However, these formulae are rather cumbersome and time consum- 
ing in use. For interpolation with more than one variable, formulae in terms of (1 
the tabular entries f;; directly instead of differences are much simpler to use. 
SALZER! has derived such a formula for double-forward interpolation in terms of 
fiz. This formula works satisfactorily for interpolation near the head of a table. w 


For the purpose of interpolating at other parts of a table, double interpolation 
formulae for other than double-forward interpolation formula have been derived. 
Expressions for partial derivatives of different orders in terms of both double 
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differences and f;; directly, including those needed in the numerical solution of 
partial differential equations, are also worked out. 

2. Double Interpolation Formulae. For making a double interpolation, the 
tabular values are usually given in a square or rectangular region. Let us divide 
this region into nine smaller regions as indicated in Fig. 1 by Roman numerals. 
In order to interpolate values in all nine regions, nine different interpolation 


II VII IV y 
VIII Vv IX 
I VI Ill 


Fic. 1. A tabular values block 


formulae are needed. These nine formulae can be further classified into five types, 
i.e., double-forward, semi-forward, double-backward, semi-central, and double- 
central interpolation formulae, as shown in Table 1. The double interpolation 
formulae now available are mostly given as functions of double differences. The 
interpolation formula in terms of double forward differences is known as Bier- 
mann’s formula or the double Gregory-Newton formula. There are three double 


TABLE I. Different types of Interpolation Formulae 


Interpolation Formula Type 


Forward x, Backward y ; 

Backward x, Forward y Semi-Forward 
Backward x, Backward y Double Backward 
Central Differences Double Central 
Central x, Forward y 
Central x, Backward y 
Forward x, Central y 
Backward x, Central y 


Semi-Central 


central difference formulae: Stirling-Stirling, Bessel-Bessel and Stirling-Bessel. 
The other types of the double interpolation formulae as functions ef mixed double 
differences may also be derived easily. 

However, these formulae in terms of double differences are rather cumbersome 
and time consuming in use. For interpolation with more than one variable, for- 
mulae in terms of the tabular entries f;; directly instead of double differences 
should be much simpler to use. Salzer'* has derived such a formula for double 
forward interpolation in terms of f;; directly as given in the following equation 


where = =—*", = and h and are the tabular intervals. 


denotes with (5) 
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This formula works satisfactorily for interpolation near the head of a table or in 
region 1. For the purpose of interpolation in other regions, similar double inter- 
polation formulae of f;; for other than double forward interpolation formula 
are required. 

It is found that by correctly changing the signs of u and v as well as the 
corresponding subscripts of f;; in Salzer’s formula, the interpolation formulae 
which apply in regions II, III, and IV can be obtained. For example, in order to 
get a semi-forward interpolation formula which applies in region II, i.e., forward 
x and backward y, one simply substitutes v into —v and f,; into f;,_; in equation 
(1), and obtains 


Similarly, the other semi-forward interpolation formula which applies in region 
III can be obtained by substituting —u by u, and fi; by f_;,; in equation (1). If 
one changes both u to —u and v to —v and at the same time fj; to f_; —;, a double- 
backward interpolation formula results. 

These semi-forward and double backward interpolation formulae of f;; thus 
obtained may be proved to give the same results as those obtained from the corre- 
sponding interpolation formulae in terms of double differences. 

Any semi-central interpolation formula, for example central x and forward y, 
should make the full use of the tabular entries of both f;; and f_; ;. One of the 
simplest ways of doing this is the following: 


v) = 3Lf(u, 2) + f(—u, »)]. 


This semi-central interpolation formula can be proved to be correct by the same 
procedure as employed by Salzer! in proving the double forward interpolation 
formula, provided one can change f(u,v) into f(u, —v) whenever v < 0, and 
f(u, v) into f(—u, v) whenever u < 0. The other three semi-central interpolation 
formulae can be obtained similarly. Finally, the double-central interpolation 
formula may be written as follows: 


f(u, +) = + f(u, —2) + f(—u, 2) + f(—u, —2)]. 


Even though the interpolations made by using these semi-central and double- 
central interpolation formulae are very close to those calculated from the double 
interpolation formulae involving Stirling and Bessel central differences, the results 
are not exactly the same. Since the interpolation formulae of f;; use the tabular 
entries more symmetrically than those formulae of central differences do, it may 
be expected that the interpolation formulae of f;; give more convergent results. 

3. Double Differentiation Formulae. In solving a differential equation nu- 
merically, either by an iterative method or by the step-ahead method, all the 
derivative terms of the differential equation should be expressed in terms of finite 
differences. Formulae for expressing the derivatives of different orders for one 
independent variable in terms of finite differences have been completely derived 
for forward, backward, and central differences.* There is thus no difficulty in 
replacing an ordinary differential equation by a difference equation. 
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However, the differential equations involved in solving most physical problems 
contain two or more independent variables. For solving partial differential equa- 
tions with two independent variables numerically, formulae for expressing partial 
derivatives by double differences are required. In order to apply to all parts of a 
field of interest, the partial derivatives should be expressed by different combina- 
tions of double differences. These expressions have been derived and the results 
are given in Tables II, III, and IV. 

In order to avoid the computations of double differences, it again may be 
worth while to work out some expressions of the partial derivatives in terms of 
tabular entries f;; directly instead of double differences. These expressions can be 
obtained from the corresponding double interpolation formulae of f;; previously 
derived. This is done by differentiating the interpolation formulae of f;; with 
respect to u and v the required number of times according to the kind and order 
of the partial derivative of interest and then letting u = 0 and v = 0. A general 
expression for the partial derivative of f(u, v) of order p + q is obtained as follows, 
where zeros are to be factored out in accordance with usual continuity convention : 


n n—i—j—1 1+b 1 j-l 
—1)i+d 
where p and g can be any integer. ‘‘’’’ on the summation sign indicates that no 
two i, terms can be identical. For instance, when i = 3 and a = 2 


The general expression of the derivative of f(u, —v) is similar to equation (2) 


except that (—1)! is removed and a (—1)™ term appears in front of (‘) and 


also fi; changes to f;_;. Similarly, the expression for the partial derivative of 
f(—u, v) can be obtained from equation (2) by removing (—1)* and placing a 


(—1)* term in front of (‘) , and at the same time changing f;; to f_; ;. If both 


of these processes are carried out, an expression of the partial derivative of 
f(—u, —v) is obtained. The partial derivative of f(-:u,v) then takes the form 
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DOUBLE INTERPOLATION 
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The expressions of the partial derivatives up to the third order have been written 
out explicitly as shown in Tables V and VI. 
Kvuo-Cau Ho 
University of Florida 
Gainesville, Florid 
1H. E. Sauzer, Jn, Math. and Physics, v. 26, 1948, Pp 294-305. 


2H. E. Sauzer, Amer. Math. Soc., Bull., v. 51, 1945, p. 279-280. 
3 Nautical Almanac, 1937, p. 802. 
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The First Published Table of Logarithms to the Base Ten 


The first two lines of the title-page of this Table. LOGARITHMORVM/ 
Chillias Prima., are followed by 26 further lines. The full title is given in JAMES 
HENDERSON, Biblioteca Tabularum Mathematicarum, 1926, p. 30; a facsimile of 
the title-page is given in A. J. THompson, Logarithmetica Britannica, part IX, 
1924; or in his completed work, frontispiece to volume 2, 1954. There is nothing 
on the title-page to indicate the author, or the date and place of publication of 
the little 16-page pamphlet, which was, apparently, privately printed. On page 
122 of Joun Warp, The Lives of the Professors of Gresham College, London, 1740, 
there is the following quotation from a letter, dated 6 December 1617, written 
by Sir Henry Bourchier to Dr. Usher: ‘Our kind friend Mr. Briggs hath lately 
published a supplement to the most excellent table of logarithms, which I presume 
he sent to you.’”’ Thus in connection with other facts the author was revealed to 
be Henry BriGcGs (1561-1631) and his table was printed in the latter part of 1617. 
John Napier (1550-1617) died in the previous April. Briggs had visited Napier 
in 1615, spent a month with him in 1616, and planned to show him this Table in 
the summer of 1617. 

Copies of the Table are excessively rare; the only copies known to exist are 
two in the British Museum: (a) with press mark c.54e 10(1); (b) a copy in the 
Museum’s Manuscript Room; and (c) a copy in the Savilian Library, Oxford. 
(a) and (c) are bound up with Edmund Gunter’s Canon Triangulorum, 1623, so 
that the original may have been trimmed; its present size is 9.3 X 15.5 cm. A 
photostat copy of (a) is a recent acquisition of the Library of Brown University. 

In this Table are given log N, N = [1(1)1000; 14D], with the first four num- 
bers of first differences, rounded for N = 500(1)1000. Characteristics are sepa- 
rated from the decimal parts by lines. The accuracy of this Table is very extra- 
ordinary. Every entry of the Table was compared with A. J. Thompson’s log N, 
N = [1(1)1000; 21D] and the only errors were 153 in the fourteenth decimal 
place: 150 unit errors, and 3 two-unit errors at N = 154, 239, 863. 

In contrast to this, when we turn to Briggs’ remarkable Arithmetica Logarith- 
mica, London, 1624, giving log N for N = [1(1)20 000, 90 000(1)10 000; 14D] 
with difference throughout, we find, for N = 1(1)1000, no less than 19 errors; 
18 two-units in the fourteenth decimal place and one serious error of 6 units in 
the seventh decimal place. Of the 150 unit errors in the 1617 publication the same 
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unit errors occur in 121 of the entries of the 1624 volume. In neither table did we 
make any checking of difference entries. 

The idea of constructing a table in which the logarithm of unity was zero 
originated with Napier. Napier and Briggs never thought of logarithms as ex- 
ponents of a base. An excellent exposition of their ideas is given in G. A. Gipson, 
“‘Napier’s logarithms and the change to Briggs’s logarithms,” p. 111-137 of C. G. 
Knott, Napier Tercentenary Memorial Volume, London, 1915; see also, H. S. 
CarsLaw, ‘‘The discovery of logarithms by Napier,” Mathematical Gazette, v. 8, 
1915, p. 76-84, 115-119. It was not till considerably later that our modern defini- 
tion of a logarithm as an exponent was put forward by such mathematicians as 
Davip GREGORY, 1684; Wm. GARDINER, 1742; LEONARD Ever, 1748, 1770. 

R. C. ARCHIBALD 
Brown University 
Providence, R. I. 


On the Numerical Integration of Functions Tabulated 
in Logarithmic Form* 

In a number of physical problems which can be described by differential 
equations, it occurs that one or more of the variables show a range of several 
orders of magnitude. In such cases, it is convenient to tabulate the variables in 
logarithmic form. The purpose of this note is to show that the usual finite-differ- 
ence formulae for numerical integration can be easily adapted to integrate a 
function when only its logarithm is given. 

We shall limit our attention to the Newton-Gregory (backward-difference) 
formula, which is the most commonly used in the hand-integration of differential 
equations. Let y’ be the function to be integrated with respect to the independent 
variable x between the limits x» — h and xo, where h is the interval of tabulation. 
Writing x = xo + hm, we can approximate In y’ by the “Newton-backward” 
interpolating polynomial f(m) as follows: 


1) Ing = som) = fot 


= fot 4%, + truncation error. 
j=l 
G= 1, 2, 3, m). 


We now want to obtain an integration formula of the type 


(2) in(t yas) =In f "exp f(m)dm 


=fo+in {r+ 0%, 


n ges 2 
= fot + N"A_, + P 
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Dropping the subscripts under the A®’s, which are always — $j, and making n 
successively equal to 1, 2, etc., we have 
1 
exp (N’A’) = fie (mA’')dm = (1 — e*’) = F, 
(3) 
0 
exp (N’A’ + N”A”) = free [mA’ + 4(m? + m)A’" = F", etc. 


The individual coefficients N can be computed from 
(4) N@A® = In F® — In FH”, 


The evaluation of the coefficients NV in a power series of A is quite laborious 
for 7 > 2; it can be easily shown, however, that the first term of the series is 
always the corresponding Newton-Gregory coefficient : 


1 1 


1 
” 2+ 244 2880 
12 1440 720 181440 


1 
N mt ‘ 


While the presence of the large term + yA’ in the expression for N’ makes 
it necessary to tabulate this coefficient (or, rather, the product N’A’) for practical 
work, it appears that the first term (i.e., the Newton-Gregory coefficient) is a 
sufficient approximation of N“ for 7 > 2. The series for N’’ shows that even in 
the extremely unlikely case of A’ = A” = 1, we have N” = —1/12.3, and nu- 
merical integration of Equation (3) shows that when A’, A”, and A’”’ vary be- 
tween 0 and 1, N’” varies between —1/24.0 and —1/24.5. In practice, the interval 
of tabulation is always chosen such that |A’| <1 (A’ = 1 corresponds to an 
increment by a factor of e in the variable from one tabular value to the next!), 
|A”’| < 1/10, |A’”’| <|A”|, etc. Under such circumstances, the coefficients N’”’ 
and N’” do not differ from the Newton-Gregory coefficients —1/12 and —1/24 
by more than 1% and it stands to reason (although convergence proofs have 
not been undertaken) that the same will be true of the coefficients of the higher 
differences. 

The method here described has been used with great success in the integration 
of the differential equations governing the motion and the loss of mass of meteors. 
Decimal, rather than natural, logarithms have been used throughout, and a 
small table of N’A’ to four places, computed for the range from —0.40 to +0.40 
in A’ at 0.01-intervals has proved more than adequate. The table is not reproduced 
here, since it can be duplicated by any computer with very little effort. 

The advantage of using the above formula should be quite obvious in the case 
of near-exponential functions, for which the step of integration can be made 
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much larger, without any increase in the truncation error, when the logarithm 
rather than the natural value of the function is tabulated. 

Luici JACCHIA 
Harvard College Observatory 
Cambridge, Massachusetts 


* Work performed under Contract DA 19-020-ORD-2556 with the Office of Ordnance Re- 
search, U. S. Army. 


Note on a Logarithm Algorithm 


In a recent note,! D. SHANKs developed a well-known algorithm for the compu- 
tation of logarithms,” in a way particularly suitable for use on automatic com- 
puting machines. In what follows I should like to point out that, if we use the 
value of » = In do (the natural logarithm of ao), the number of operations neces- 
sary for this computation can be cut down considerably in replacing about a 
third of the single steps indicated by Shanks by one division and one addition. 

We assume as in the paper quoted that a) > a; > 1. Tocompute A = log,, a1, 
we determine a sequence of numbers dz, 3, d4, ---, (@; > 1) and a sequence of 
positive integers m1, m2, ---, by the relations 


We then have 
1 1 


if we stop at the calculation of m; we have an approximate value of \ by taking 
the ith convergent: 


& 


We will show that we have for 7; the formula 


— 1 
(I) n=n* +p, = 
where the error term p; can be estimated by 
(II) wOm? < < 
as soon as we have 
(IIT) | u/Qi < 1.7933 ---. 


For ao = 10 we have w = pio + 2.3026 and (III) is certainly satisfied from 
4 = 2 on. For ao = 2 we have 


B= pe = 0.6931. 
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Thus if we want, e.g., to compute log: a; with an error less than 10~°? it is sufficient 
to stop with the computation of n; as soon as Q; exceeds 10*. In determining n; 
by successive division we will already have the value of a; = a;_;/a;"‘ and then 
apply 

In order to prove (I)—(III), observe that we have for all positive integers 


2-1 
ao 


Indeed, since P; = 1,Q:1 = 1, P2 = m2,Q2 = mymz + 1, our formulae are verified 
immediately for 
ao a, 
a3; = 
a," a," ao" 


From there on we proceed by induction, using the well-known relations 
P; = + Pi-e, Qi = niQi-1 + Qi-2. 
If we write now e = (—1)‘ we have from (1), for an even i = 2p: 


= CQ = 


that is, 
(2) = 


this formula also holds in the case of an odd 7, where we have « = — 1. 
From (2) we have now 


(3) — 1 = + 
where @ is the value of the series 


Therefore we have certainly |@|< 1, if |y| does not exceed the positive root yo 
of the equation 


v=1+p+P'; 
we find 
vo = 1.7933. 
Assuming now that we have 
(4) < vo, 


we obtain (I), with p; = — 0(uQm?), on dividing (3) by euQ;. The estimates (II) 
follow then immediately as we have |n;:| < 1/Q?, |@| < 1. On the other hand, if 
(III) is satisfied, so is (4), since Q;|n;| < 1/Q;. 

We give for a) = 10, a; = 2 the exact values of the ; and the values 7;* 
computed by (I), in the last two columns of the following table for 7 = 1, ---, 7. 
The three first columns of this table contain the corresponding values of Q;, a; 
and y;. Although the values of the a; are given in the same example treated by 
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Shanks with 9 decimals fori = 1, ---, 6, we have recomputed them and computed 
further a; and as on the 10 X 10 Friden Desk Calculator without double precision 
or any similar artifice. 

In computing P;/Q; + »;* it is better not to compute the second term sepa- 
rately from (1) but to use the formula 


+ = [>. + (- 1) su —t1/o, 
in replacing in this way two divisions by one only. 
a; Mi ni* 

.33333 33333 333 — 3.619 —3.2 1 
10 1.25 .30000 00000 000 .042 x 10° 1.03 1 
93 .30107 52688 1720 — 4.549 — 4.527 1 

1 

1 


2136 1.00104 1546 99625 468 


3 x 
i x 
x 
00974 1958  .30102 04081 6327 10-* 95875 xX 

x 1 
3 x 1 
13301 1.00016 2900  .30102 99977 4453 —2 x 1 


.00433 6278  .30103 09278 3505 —9. 


We see from this table that, while the error of us is 3.3 X 10-*, the improved 
method gives already us + 4* with the error 2.1 X 10-*. On the other hand in 
using all computations which give us we obtain at once a; and from there on 
obtain us + ye* by (I); this is in error by about 3.2 X 10-”. If we go on toi = 7 
the error of yu; is 2.2 X 10-*, while since we obtain at once as, we get from (I) 
the approximation y; + 77* the error of which is about 9.0 x 10-*. Here, how- 
ever, the result is already influenced by the rounding off error in our value of as. 
If we use instead the value of as obtained by double precision : 1.00006 37223 565, 
we obtain for ;* the improved value 4;* = — 2.0806156 X 10~ and the error of 
uz + H7* is then 6.6 X 10-“. We see that it is hardly worth-while to go over to 
double precision unless we need essentially more than 10 decimals. 

As the essential steps in the original method consist in the successive deter- 
mination of integers, it is clear that this method is not very sensitive with respect 
to the rounding off. It appears a little unexpected that the same is true to a very 
great extent for the improvement of this method proposed here. It turns out that 
the denominators 4Q; in the formula (I) are sufficiently great to counterbalance 
to a certain extent the rounding off errors in the a;. 

It may be finally remarked that from a certain point of view the conclusions 
of D. Shanks concerning the linear character of convergence are perhaps slightly 
too optimistic. As a matter of fact the amount of work will not be measured by 
the number of ‘‘cycles,”’ i.e., the number of partial quotients of the continued 
fraction used, but the total number of divisions; and this will be measured not by 
the number of the »; but by their sum m; + m_. + --- + ;. About this sum, 
however, KHINTCHINE has proved‘ that its order “in the average”’ is not that of 
4 but that of zlg 7. The “linearity of convergence’’ will probably fall off in the 
long run in the “general case.’’® 

A. M. Ostrowsk1 


Mathematical Institute of the University 
Basle, Switzerland 
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1 DANIEL SHanks, “‘A Logarithm Algorithm,” MTAC, v. 8, 1954, p. 60-64. 

2 See, e.g., J. TROPFKE, Geschichte der Elementar- Mathematik. Bd. II: 3. Aufi., 1933, p. 241-242. 

In applying (I) we can of course replace the division by yu by the multiplication with 
M = 1/p = logay é, i.e., with the module of the logarithms to the base do. In this way we replace 
two divisions indicated in (I) by one division and one multiplication. 
pe A. . Metrische Kettenbruchprobleme. Comp. Math., v. 1, 1935, p. 360-382, espe- 

ially p. 376, 4 

$ Br. D. SHANKS to whom this note was submitted in manuscript made the very interesting 

remark that the use of » = Ina can be eliminated from my formula (I) in replacing u by 
Pi(a; — 1) + — 1) + 


However, in this case, the exact error estimate in (II) has to be slightly changed. 


Iterative Procedures for Taking Roots Based on Square Roots 


Introduction. Several years ago, on the Model 1 C.P.C. at the Los Alamos 
Scientific Laboratory, there was incorporated by Dr. R. H. Stark a routine 
for taking square roots which proceeded at the same rate as the basic card 
feed operation. In connection with such an operation it is natural to ask if cer- 
tain basic calculations might not be simplified by adding the square root as an 
operation. Today the existence of the Friden desk calculator, which also takes 
square roots directly, suggests that the question continues to have some cogency. 

In connection with the Los Alamos C.P.C., we devised a procedure for taking 
cube roots and other roots which we present here for its conceivable usefulness 
and interest. 

The Newton-Raphson iteration procedure for approximating roots of an equa- 
f(x) 
(x) 
gives an iterative procedure for finding roots starting with a guess xp. 

To take the cube root of a positive number m, then, we may consider equa- 
tions x* — n = 0, x*? — ni? = 0, x*/4 — n'/4 = 0. The corresponding respective 


iteration functions are gi(x) = 1/3 (5 + 2) » B2(x) = 1/3(2n*x-1/2 + x) and 


tion f(x) = 0 establishes a function g(x) = x — 


such that = g(xu) 


gs(x) = 1/3(4n'/4x/4 — x). Now either of the last two functions are not feasibly 
used unless one can readily take square roots or fourth roots. In general, they 
are more costly than the usual iteration g(x). However, if one has no preference 
between, say, square roots and other arithmetic operations, then either g2(x) or 
gs(x) will give an iteration which tends to converge in fewer steps than the 
standard one. For example, we have found cases in which the second routine was 
adequate in four iterations, whereas the g; routine would require five starting with 
the same guess. 

To handle negative roots one may assign the sign of to the first guess. We 
found it convenient to guess ¥ in the C.P.C. which found ¥n automatically. In 
8(nx)!8 — x 


general, one may take the seventh root of m, say, by using g(x) = . 


1/8 
and the ninth root by using g(x) = 1/9 (s (2) + ) , and so on. 


Preston C. HAMMER 
University of Wisconsin 
Madison, Wisconsin 
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A Note on the Electronic Computer at Rothamsted 


The “401” experimental electronic computer built by Elliott Brothers under 
contract from the National Research Development Corporation is now housed 
in the Statistical Department at Rothamsted, machine time being shared on a 
50-50 basis between the N.R.D.C. and the Department. The Department is 
applying the computer to probiems which arise in statistical research and its 
applications to agriculture, but time on the machine is available through the 
N.R.D.C. for general computations. The use of the computer specifically for the 
former work is a new venture. Prior to Rothamsted’s acquiring the 401 little 
was known of the advantages and difficulties of applying electronic computers 
to statistics. Consequently, before being able to take full advantage of this power- 
ful aid to computa‘ion, it was necessary to obtain experience not only of the 
appropriate programming techniques but also of the methods of adapting the 
rather special numerical techniques. By using a prototype computer the relevant 
experience is being gained and at the same time the characteristics that a perma- 
nent machine in this field ought to have are becoming apparent. This note gives 
a short description of the machine and an indication of the type of work being 
carried out. 

The 401 is a serial computer with a word length of 32 binary digits and a word 
time of 100us. The main store is a magnetic disc of 23 tracks, each holding 128 
words; 7 of the tracks are always immediately available to the computer, but 
only one of the remaining 16 is accessible at any given time, switching between 
these tracks being by means of high-speed relays. There are 5 single word im- 
mediate-access registers; one is the accumulator which can be coupled with the 
second register for double length working ; the other 3 registers can each be used 
to modify orders as well as for temporary storage of numbers. Input is by 5-hole 
punched tape and output by either electrical typewriter or teleprinter-punch ; 
punch-card input is shortly to be incorporated. The machine has a two-address 
code to allow for optimum programming. Operations possible include addition, 
subtraction, multiplication with and without round off, collation, non-equivalence, 
discrimination on zero or negative values, and left and right shifting of single or 
double length numbers. Addition and subtraction take one word time (100us) 
to carry out, a shift of m places takes m word times and multiplication 32 word 
times (3.2 m.sec.). 

A library of sub-routines is being built up and many of the standard functions 
and operations are now available. Routines for analysing randomised blocks and 
Latin square experiments have been used to obtain results for several large series 
of experiments. Work on routines for the analysis of factorial experiments is under 
way—these will allow a fast and accurate assessment of a wide variety of experi- 
ments. A programme has been developed which calculates the means, standard 
deviations, variances, covatiances and correlation coefficients for samples of 
multivariate data. Specific problems tackled include a sampling investigation of 
the “combination of probabilities” test of significance: this necessitated a very 
large amount of computation and has been successfully concluded. General 
routines developed include a ‘‘programming’’ routine and a set of routines for 
floating binary operations. The former inserts the appropriate addresses in a 
skeleton programme and punches out the final routine for immediate use, thus 
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saving a great deal of programmer’s time. The latter allows for the full 9 signifi- 
cant decimal figure accuracy and can accommodate numbers in the range 
108000, 

The encouraging start made at Rothamsted seems to indicate that electronic 
computers are well suited for dealing with statistical problems, but undoubtedly 
new numerical and programming techniques are necessary to meet the special 
requirements of such work. 


S. Lreron 
Rothamsted Experimental Station 
Harpenden 
Herts., England 


REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


This section contains reviews and descriptions of tables and of publications 
which in some explicit way have general interest in connection with tables or 
computation. Those noted as ‘‘Deposited in the UMT File’’ have been deposited 
in the Unpublished Mathematical Tables file maintained by the Chairman of 
the Editorial Committee; they are available there for reference. In many cases 
the author of these tables has a limited number of copies available for distribution. 
Authors of works containing mathematical tables of general interest or which 
otherwise fall into the classes noted or reviewed here are urged to submit a copy 
to the Chairman of the Editorial Committee. 


31[A, F].—RSMTC Table of Binomial Coefficients. Royal Soc. Math. Tables, 
v. 3. Edited by J. C. P. MiLLer. Cambridge 1954; viii + 162 p. 21.5 & 28.0 
cm. Price $6.50. 


This table gives exact values of binomial coefficients 


for positive integer values of k and m. The main table (p. 2-103) is for m < 200. 


Because 
2) 
k n—k 
it is unnecessary to print values for k > m/2. Actually k is taken less than or 
equal to (m + 1)/2 plus whatever additional values of & fill out the line. The 
central binomial coefficient in case m is even, or the two central coefficients in 
case m is odd, are printed in bold face type. 
The rest of the table is for  < 5000 but for k < L limited as follows 


n 
200-500 12 
500-1000 11 

1000-2000 5 


2000-5000 3 


a 
T 
d 
h 


42 
T 
| 
I 
: 
‘ 


REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 71 


Thus the table contains in particular the first 5000 triangular and tetrahedral 
numbers. 

The arrangement of the table is such that the coefficients ( 7 ) and ( k r4 n 
are in general adjacent entries in the same line; that is, column headings are k. 


The opposite arrangement in which 4) (k = 0,1, 2,---) form a column of 


data would be perhaps more convenient for most purposes but would, no doubt, 
have rendered a difficult printing problem more difficult. On page 1 there is an 
index by means of which the reader can find any particular binomial coefficient. 

The printing is by Cambridge University Press and is superlative, as always. 

This table is intended primarily for number theorists whose interests in bi- 
nomial coefficients are such that exact values are required. For the more practical 
man who wants no more than a half dozen significant digits the table will no 
doubt prove irritating in spite of the fact that the printer has supplied auxiliary 
figures to help him truncate the exact values and to express them in floating 
decimal notation. However the practical man may also have his curiosity aroused 
by inspection of the table. He may observe for example that for a fixed m the 
number of odd binomial coefficients is always a power of 2. For such reasons the 
reviewer would recommend the table for general use. These are the most extensive 
tables of this function yet published. They are the result of many collaborating 
calculators. 

D. H. LEHMER 

Univ. of California 
Berkeley, Calif. 


32[B ].—ApmrrALty RESEARCH LABorATORY, “Tables of » = {m? + ¥— 1},” 
A.R.L./T.7/Maths 2.7, April 1954, 4 p., Teddington, Middlesex, England. 


The tabulation to six decimal places with second differences of |y|, arg 2, 
R(n) and I(n) is given for m = 0(.02)3(.1) up to the point where » = m + 4m-', 
arg 7 = 4m~*, R(n) = m + $m, I(n) = 3m, to the accuracy of the table. 

I. A. STEGUN 
NBSCL 


33[D].—NBS Tables of Circular and Hyperbolic Sines and Cosines for Radian 
Arguments. Appl. Math. Ser. No. 36. U. S. Gov. Printing Office, Washington, 
1953. x + 407 p. Price $3.00. 


This is the third edition, published in 1953, of this set of tables originally 
published in 1939. See RMT 89[D, E], MTAC, v. 1, 1943-45, p. 45, and Note 
114, v. 4, 1950, p. 123. The present edition corrects two inconsistencies of format 
and six printing errors of the first edition ; five of these printing errors have been 
corrected in the second edition. 

The only other reported change is an extension of a supplementary table, 
which now expresses degrees, minutes, and seconds, in terms of radians to 10 
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decimal places and selected values of radians in terms of degrees, minutes, and 
seconds, to an accuracy of 0.000005 second. 


34(D].—L. W. Potiak, Eight-Place Supplement to Harmonic Analysis and Syn- 
thesis Schedules for Three to One Hundred Equidistant Values of Empiric Func- 
tions. Second Edition, Dublin Institute for Advanced Studies, School of 
Cosmic Physics, Geophysical Memoirs No. 1, Dublin, 1954. Price £1 2s. 6d. 


This is a second edition to the work issued in 1949, RMT 847[K ] reviewed in 
MTAC, v. 5, 1951, p. 19-21. The work was written to be associated with “‘All 
term guide for harmonic analysis and synthesis,” see RMT 899[K ], MTAC, v. 5, 
1951, p. 149. However, directions are included for its use independently of the 
earlier “‘Guide.” 

In the present issue the accuracy of angles given in degrves, minutes, seconds, 
and fractions of a second, columns headed 7z and dz’ has been carried to ten decimal 
places of a second. Some checks of accuracy of the earlier work are reported, and 
it is reported no mistakes were discovered. 

TF. 


35[E].—AECD-3497, S. FRANKEL and E. Netson, ‘Methods of treatment of 
displacement integral equations,’’ Sept. 1953. Available from Office of Tech- 
nical Services, Department of Commerce, Washington 25, D. C. 


This contains on p. 78-82, a table of k = k(f) defined by 


1 arctanh k 
= < < = = 
i<f<90, i+f 
1 arctan k 
<7 = 
0<f<s, i+f 


f = —1(.01)2(.1)5(.5)10. The table is to 6D, with first differences given. 
A polynomial of degree 5 in (1 + f)~' which gives 4D in the range 8 < f < ~, 
is given. 


36[K].—Lesiiz E. Suwon and Frank E. Grusss, Tables of the Cumulative 
Binomial Probabilities. Ballistic Research Laboratories, Ordnance Corps. 
Pamphlet ORDP20-1, Aberdeen Proving Ground, 1952. 577 p. + viii, 12 x 9 
inches, photo-offset from typescript. Available at U. S. Gov. Printing Office, 
Washington, D. C. 


This table lists 


Plan p) = py 


for c = 1(1) --- to P < 10-7, m = 1(1)150, and p = .01(.01)50.00. This carries 
nm much farther than other available tables. 
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The introduction and foreword explain the use of the tables and give three 
illustrative examples. It is pointed out that the Incomplete Beta function 
I,(c,m — ¢ +1) = P(c, n, p) can be read from the table. 

The type is small but is generously spaced: each value of m has an integral 
number of pages for its own. 


H. CAMPAIGNE 
National Security Agency 
Washington, D. C. 


37[L].—L. Fox, A Short Table for Bessel Functions of Integer Order and Large 
Arguments. Cambridge, 1954. Royal Society Shorter Mathematical Tables, 
No. 3. 28 p. Price $1.25. 27.8 X 21.7 cm. 


The present tables are intended to supplement the tables of Bessel functions 
of integral order, J,, Y,, J, and K, previously prepared by the British Association 
Mathematical Tables Committee. In the earlier volumes the functions were 
tabulated for m varying from 0 to 20 and x in the range from 0 to 20 or 25. Here 
the range of x extends from 20 to «. However, the tabular argument is 1/x so 
that use of the tables requires calculation of a reciprocal before entering. This 
device permits coverage of the range by means of fifty tabular values, namely 
1/x = z = 0(.001).05. In the case of J,(x), Y,(x) the usual auxiliary functions 
P(x), Q(x), for which 


Jn(x) ~ [P.(x) cos (x — 4nx — — Q,(x) sin (x — 4nx — }x)], 


Y,.(x) [P.(x) sin (x — 4nx — 1/4x) + Q,(x) cos (x — — }x)], 


are tabulated. In the case of I,,(x), K,(x), the auxiliary functions F,(x), Y,(x) 
such that 


ez 


F,(x), K,a(x) ~ vz e~*G,,(x) 


are tabulated. P, and Q, are given to nine decimals for m = 0(1)9 and to eight 
for n > 10. F, and G, are given to nine decimals for n = 0(1)9 while In F,, and 
In G, are given to eight decimals for m > 10. Modified second and fourth central 
differences are given for interpolation. 

P,, and Q, were computed for m = 0 and m = 1 for z = 0(.005).07 and the 
recurrence relations 


I,(x) ~ 


Pay — Pai = —2nzQ,, 


were used to generate the values for those arguments for m = 2(1)10. The values 
were then subtabulated at intervals of .001 for m = 0(1)10. For n > 9 subtabula- 
tion was not always convenient and the values at intervals of .001 were generated 
by recurrence. Similar methods were used for the computation of F, and G,. 
The tables were reproduced by photo-offset from a manuscript prepared on a 
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card-operated typewriter. From the numerous checks applied in the course of 
preparation one may feel confident that these tables are free from error. 


MILTON ABRAMOWITZ 
National Bureau of Standards 
Washington, D. C. 


British Association Mathematical Tables, vol. 6. Bessel Functions, Part I: Functions of orders 
Zero University Press, 1937, reprinted 1950. [See RMT 179[L], MTAC, 
v. 1, 1944, p. 361. 

British Association Mathematical Tables, vol. 10. Bessel Functions, Part II: Functions of 
port integer order. Cambridge University Press, 1952. [See RMT 1087[L], MTAC, v. 7, 1953, 
p. 97. 


38[L].—NBSCL, Table of the Gamma Function for Complex Arguments. NBS 
Applied Mathematics Series No. 34. U. S. Gov. Printing Office, Washington, 
D. C., 1954, xvi + 105 p., 22 * 20 cm. Price $2.00. 


Values of the gamma function for complex arguments are required for work 
in atomic and nuclear physics, engineering, and elsewhere, and the publication 
of a reliable basic table of this function is to be welcomed. The present volume 
gives 12D values of the real and imaginary parts of log, T(x + zy) forx = 0(.1)10, 
y = 0(.1)10. The values for negative y may be obtained simply by changing 
the sign of the imaginary part, the values for negative x from the relation 
T(z)T'(1 — 2) = wcosec rz. To facilitate computation from this last formula, 
auxiliary tables of sin rx, cos rx, sinh rx, cosh rx to 15D or 15S for x = 0(.1)10 
are included. Thus, the tables given in this volume may be said to cover the 
region — 9 < x < 10, — 10 < y < 10. The effective range of the tables may be 
doubled in each direction by an application of the duplication formula; and out- 
side the so increased range, the asymptotic expansion of log '(z) may be used. 

Key values of log, ' (x + iy) for x = 9(.1)10, y =0(.1)10 were computed to 
15D from Stirling’s expansion, the values for 0 < x < 8.9 were obtained from the 
key values by application of the recurrence relation ['(z) = z"IT'(z + 1). A de- 
tailed description of the computation is given in the Introduction by H. E. 
SALZER which contains also notes on the purpose and scope of the table, on some 
properties of the gamma function, and on direct and inverse interpolation in the 
table. Also included is a bibliography listing 11 other tables of gamma functions 
for complex arguments, 16 auxiliary tables, and 16 references to books and papers 
containing relevant material on the gamma function. 

The technical staff of the Mathematical Tables Project at the start of this 
undertaking consisted of A. N. Lowan, Chief, Mitton ABRAMOWITZ, GERTRUDE 
BLAaNcH, ABRAHAM HILLMAN, WILLIAM HORENSTEIN, MEYER KARLIN, JACK 
LADERMAN, IDA Ruopes, H. E. SAuzer, and IRENE SteGuN. The key values for 
9 < x < 10 were computed by hand on desk calculators by RutH Capuano, 
LEONA FREEMAN, ABRAHAM GROSSMAN, and RuTH ZUCKER. The extension to 
0 < x < 8.9 was done largely on punch-card machinery, by MILTon STEIn, who 
was also responsible for preparing and checking the final manuscript. 

A. E. 


39[L ].—AKapemrrA Nauk SSSR. [Institute of Exact Mechanics and Computa- 
tional Techniques. ] Tablitsy Integralnogo sinusa i kosinusa [Tables of the sine 
and cosine integral |. Moscow, 1954. 473 p. 25.7 X 19.7 cm. 43.75 rubles. 


1 
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These tables give, to 7D, 
Ci(x) = f * 1 cos tdt 


for x = 0(.0001)2(.001)10(.005)100. 

There is a one page table of C;(x) — In x for x = 0(.0001).0099. These are 
rounded values taken from the 1940, 1942 NYMTP tables of these functions,' 
except that the range 10 < x < 100, with the old interval of .01, has been sub- 
tabulated to .005. This was apparently done by hand, punched on cards, and 
checked with a tabulator. In most places no differences are given since linear 
interpolation suffices. When they are not negligible, second differences are given. 
A one page table of $#(1 — #) and a nomogram for interpolating with A* are 
supplied as inserts. The table is well printed on good quality paper. With good 
proof reading this could be a very accurate table. If one does not need the ex- 
ponential integral, or more than 7D accuracy, this volume is a handy replacement 
for the three NYMTP volumes. 


D. H. LEnMER 
Univ. of California 
Berkeley, Calif. 


1NYMTP. Tables of Sine, Cosine and Exponential Integrals. V. 1, 2. New York, 1940. Table 
of Sine and Cosine Integrals. New York, 1942. 


40[L].—AKapemiiA Nauk SSSR. [Institute of Exact Mechanics and Computa- 
tional Techniques. ] Tabliisy Integralov Frenelya [Tables of Fresnel Integrals}. 
Moscow, 1953. 271 p. 16.2 X 25.8 cm. 23.50 rubles. 


These fundamentally new tables of the Fresnel integrals 


S(x) = f ein f 


C(x) = f cos = ; f 


give 7D values of these functions for 
x = 0(.001)25 
together with average second differences 
${A*f(x — Ax) + A*f(x)}. 


Two short tables near x = 0 give 7S values (rather than the 7D values of the 
main table). These are tables of 


S(x) for x = 0(.001).58 
and 


C(x) for x = 0(.001).101. 
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For interpolation there are two tables of #(1 — #)/2 for ¢ = 0(.001).5 and a nomo- 
gram for finding $¢(1 — #)A*. 

The table was found by successive quadrature performed by Simpson’s rule 
between the points x = m and x = n + 1, m = 0(1)24. These 26 key values were 
computed from the power series and asymptotic formulas for S and C. 

The most extensive previous tables of the Fresnel integrals are those of 
van WIJNGAARDEN & SCHEEN [MTAC, v. 4, p. 155] which appeared in 1949. 
These tables are to 5D for x = 0(.01)20. 

The tables are well printed on good quality paper and represent a very con- 
siderable number of man hours of effort. 

D. H. LEHMER 
Univ. of Calif. 
Berkeley, Calif. 


41[L ].—ApmrraLty REsEARCH “Table of | J| = | f sec Bei 
A.R.L./T8/Maths 2.7, April 1954, 13 p., Teddington, Middlesex, England. 


The function |J| is tabulated to four decimal places for w = 1(1)13. The 
intervals in the table have been so chosen that linear interpolation in the ¢ 
direction for the most part is adequate, or at most, the tabulated second differ- 
ences must be used. In the region around g = 90°, the auxiliary functions 
h(u),-7 and d are given, where |J| = h(u){1 + 2r cos Q + r*}# and @ (in degrees) 


180 
¢. 


I. A. STEGUN 
NBSCL 


42[L].—H. GELLMAN and JEAN TuckKER, “Tables of the functions Do(x) and 
D,(x).” Atomic Energy of Canada Limited, Chalk River, Ontario, 1954. 
51 p., 21.5 & 27.5 cm. $1.50. 


The functions tabulated in this report are 
1 (7? 
Do(x) = if d0 exp (— x cos @) 


and D,(x) = — Do'(x). They can be expressed in terms of modified Bessel func- 
tions and modified Struve functions as 


Do(x) = Lo(x), D,(x) = - I;(x) Li (x). 


Do(x) and D,(x) have been calculated on FERUT, the Ferranti Universal Digital 
Computor at the University of Toronto, and 7D values of them are given here 
for x = 0(.001)4.999. 


A. E. 
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43[L].—M. Roruman, “Tables of the integrals and differential coefficients of 
Gi(+x) and Hi(—x).” Quart. Jn. Mech. Appl. Math., v. 7, 1954. p. 379-384. 


The functions 
us 
were introduced by Scorer.' The present paper contains 7D tables of 


for x = 0(.1)10, with 7,,”. The tables were obtained by numerical integration and 
differentiation of 10D tables of Gi(x) and Hi(—<), and the results were compared 
with those obtained from the asymptotic expansions for large x. The author has 
computed values up to x = 20, which are available if required. 

A. E. 


1R. S. Scorer, ‘Numerical evaluation of integrals of the form J = and the 
tabulation of the function Gi(s) = : JF sin (us + 4u*)du,” Quart. Jn. Mech. Appl. Math., v. 3, 


1950, p. 107-112; MTAC, v. 4, 1950, p. 215. 


44[L].—M. Roruman, Table of the integrals of Ai(+x) for x = 0.00(0.01)2.00 
and x = 0.0(0.1)10.0. i + 6 p. mimeographed. Deposited in the UMT Fine. 


Ai(x) is the Airy integral of the first kind, and the tables give 


f 


to 8D with 4,,?. The integrals were computed at intervals of 0.1 using the method 
of reduced derivatives. J. C. P. MILLER’s unpublished table of reduced derivatives 
of Ai(x) was used. A shortened version of this table has been published (see RMT 
1231, MTAC, v. 8, 1954, p. 162). 

M. ROTHMAN 


Northern Polytechnic 
Holloway, London N.W. 7 


England 


45[L].—M. Rotuman, Table of Bi(+x) for x = 0.00(0.01)2.00 and Bi(—<x) for 
x = 0.00(0.01)10.00. i + 6 p. mimeographed. Deposited in :x1e UMT Fixe. 


Bi(x) is the Airy integral of the second kind, and the tables give 8D values of 
Bi(+x) with 4,,*. The computation was carried out using the method of reduced 
derivatives and 12D were carried throughout the calculations. J. C. P. MILLER’s 
manuscript tables of Bi(x) and its reduced derivatives were utilized. 

M. ROTHMAN 
Northern Polytechnic 
Holloway, London N. 7 
England 
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46[L].—M. Roruman, Table of the integrals of Bi(+x) for x = 0.00(0.01)2.00 
and of Bi(—x) for x = 0.0(0.1)10.0. i + 6 p. mimeographed. Deposited in the 


UMT 
f 
0 


8D tables, with 6,,?, of 
Bi(x) being the Airy integral of the second kind. The integrals were computed at 
intervals of 0.1 using the method of reduced derivatives. J. C. P. MILLER’s manu- 
script tables of reduced derivatives of Bi(x) were used. The integrals were also 
computed from their ascending power series for |x| < 2. 

A shortened version of this table has been published (RMT 1231, MTAC, 
v. 8, 1954, p. 162). 


M. ROTHMAN 
Northern Polytechnic 
Holloway, London N. 7 
England 


47[L ].—M. Roruman, Table of Gi(x} and its derivative for x = 0.0(0.1)25(1)75. 
i + 6 p. mimeographed. Deposited in the UMT Fie. 


A shortened version of the table of Gi’(x) has been published (see review 41 
above). The mimeographed pages give 8D values of Gi(x) and Gi’(x) with 4,,? 
(and for x < 2.5 also y*). The table of Gi’(x) was obtained by numerical differ- 
entiation of a 10D table of Gi(x). 


M. RoTHMAN 
Northern Polytechnic 
Holloway, London N. 7 
England 


48[L ].—M. Rotuman, Table of Hi(—x) and its derivative for x = 0.0(0.1)25(1)75. 
i + 6 p. mimeographed. Deposited in the UMT Fine. 
A shortened version of the table for Hi’(—x) has been published (see review 
41 above). The mimeographed pages give 8D values of Hi(—x) and Hi’(—x) with 
bn? (and for x < 2 also y*). The table of Hi’(—x) was obtained by numerical 
differentiation of a 10D table of Hi(—x). 


M. RoTHMAN 
Northern Polytechnic, 
Holloway, London N. 7 


England 


49[L ].—M. Rotuman, Tables of the integrals of Gi(x) and Hi(—x) for x = 0.0(0.1)- 
20.0. 4 p. mimeographed. Deposited in the UMT Fizz. 


A shortened version of these tables has been published (see review 41 above). 
The mimeographed pages give 8D values, with 6,,? (and in a few places also +‘), of 


f Gi(tdt, f Hi(—2)dt. 
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The integrals were obtained by numerical integration of 10D tables of Gi(x) and 
Hi(x). They were checked from x = 8 onwards by comparison with results calcu- 
lated from the asymptotic expansions. 

M. RoTHMAN 
Northern Polytechnic 
Holloway, London N. 7 
England 


50[M ].—TERENCE BUTLER and Kari PoHLHAUSEN, Tables of definite integrals 
involving Bessel functions of the first kind. WADC Technical Report 54-420, 
Wright Air Development Center, 1954. 50 p. 21.5 K 27.5 cm. 


7- is the rth positive zero of Jo(z). All values are tabulated to 5D. 

Table I. y, and J:(y,-), r = 1(1)10. This table is taken from H. T. Davis and 
W. J. Krrxuam, Bull. Amer. Math. Soc., v. 33, 1927, p. 760. 

Table II. J,(y-), 7 = 1(1)10, » = 1(1)»,, where », is the last » for which J,(y,) 
is not zero to 5D. 

Tables III-VII. 


f x? Jo™ J1"(yx)dx, 


where = 0(1)10 and r = 1(1)10 in all tables; m = 1, » = 0 in Table III; 
m = 0," = 1in Table IV; m = 2, = 0in Table V; m = 0, m = 2 in Table VI; 
m = n = 1 in Table VII. 

Tables VIII-X. 


where p = 0(1)5, r = 1(1)5, s = 1(1)5 in all tables; m = n = 0 in Table VIII; 
m = 0, n = 1 in Table IX; and m = n = 1 in Table X. 
A. E. 


51[(T, L].—Ascuer OpLer and Nevin K. Hester, Tables for Predicting the 
Performance of Fixed Bed Ion Exchange and Similar Mass Transfer Processes. 
Stanford Research Institute, Stanford, California, 1954. Multilithed 111 p., 
84 X 11 inches. Free. 


A variety of problems in the theory of non-equilibrium operation of ion- 
exchange and absorption columns in the theory of heat transfer and in the theory 
of probability may be described by the partial differential equation 


(2), is (*), = X(1 — «) — — 2). 


The solution of these equations may be written as 


» J(rs, t) 
(rs, t) + — J(s, rt)] 
1 — J(t, rs) 


ws. 
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where 


The function 1 — J(x,y) has been previously calculated by BRINKLEy, 
Epwarps, and Smita (MTAC, v. 6, 1952, p. 40). However these calculations were 
found to be insufficient for the purposes of the authors when J was near zeroor one. 

The evaluation of \ and w was carried out with an IBM 602-A using a modifi- 
cation on Onsager’s asymptotic expansion. 

Values are given to 4D and appear in three tables covering the following 
ranges. 


TableI: =0.2(0.2)1,2 s=1,2(2)8 = 0.2, 0.5, 1, 2,5 
Table II: + = 0.2(0.1)1, 1.2, 1.3, 1.4, 1.5, 2(1)5 
s = 10(5)100(10)1,000 ¢/s = 0.1(0.1)0.4(0.2)1(1)4, 6 


Table III: = 0.2(.1).9 s = 10(5)100(10)500 ¢/s chosen so that = 0.1 
and 0.9 


The report contains the derivation of the differential equations in various 
problems and a description of the method of computation as well as the tables for 
d and w. In spite of a few typographical errors (e.g., Z occurring in equation (9), 
p. 13 should be read as 2), this report is a useful collection of material for workers 
concerned with mass transfer problems. 


52[V ].—AERONAUTICAL RESEARCH CounciL, A Selection of Graphs for Use in 
Calculations of Compressible Airflow. Prepared hy the Compressible Flow 
Tables Panel (L. RosENHEAD, Chairman, W. G. BickLey, C. W. Jongs, L. 
F. Nicnotson, H. H. Pearcey, C. K. THORNHILL, R. C. ToMLtnson) of the 
ARC. Oxford, 1954. x + 115 p., 119” & 152”, 84s. net. 


This book is a companion volume to ‘‘Compressible Flow : Tables’’ (Clarendon 
Press, 1952), RMT 1093[V], v. 7, 1953, p. 103. The preface states, ‘“The object 
of both books is, briefly, to make available to engineers, physicists, and applied 
mathematicians a selection of tables and graphs likely to be of value in research 
and in calculations of the flow of air in which compressibility effects are important. 
It is hoped that the books will be useful both as an aid to design and for the 
development of new theory.” 

This handsomely printed and well thought-out book should succeed in fulfilling 
the hopes of the editors. 

There are two types of graphs in this book, corresponding to single and double 
entry tables, and for both types single-page and multiple-page graphs are given. 
The single-page graphs are intended to give a quick appreciation of the variation 
of the quantities and they can also be used for rough working. The multiple page 
graphs corresponding to single entry tables enable one to take readings which 
approach those of the tables and in some cases are much more convenient to 
use than the tables. 

In the mathematical functions plotted, y, the ratio of specific heats of air has 
been taken to be 1.4 and yu, the coefficient of viscosity for dry air at 23° C, to be 
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(483.00 + 0.25) micropoise. This coefficient is equivalent to (374 + 0.5) x 10-* 
slug/ft. sec. at the standard sea level temperature of 15° C. The upper limit of 
the Mach number of the flow has been fixed arbitrarily at 5. 

The graphs have been grouped in sections labelled A. Isentropic Flow, B. 
Normal Shocks, C. Oblique Shocks, D. Conical Flow, E. Reynolds Numbers. 
Each section contains its own introduction in which the basic equations used for 
the graphs are derived and briefly explained. The derivations do not always start 
with first principles. Nevertheless the equations plotted are clearly stated as are 
the figures in which the corresponding graphs are to be found. 

This book will be a valuable and useful tool for the engineers, physicists, and 
applied mathematicians working in compressible airflow. 

A. H. T. 


53[V].—K. G. Tapman, Tables of Flow Functions for Bodies of Revolution in 
Circular Tunnels and Jets. Armament Research Establishment Memo 25/53 
(modified). 11 mimeographed pages. 

The functions defined by the following equations occur in the computation of 
irrotational incompressible flow patterns about bodies of revolution. 


- r"P,,(s) 


o(x, —2) + 1/2 In r(1 s) = 


r—1P,,’(s) 
1/2[r(1 + — Q(x, y, », —2) = (mn + 1)! 
r*—P,,'(s) 


(n + 1)! 


where r? = x? + y*, s = x/r, P,(s) is the Legendre polynomial of order m and the 
C,,. are constants given in the memorandum to six significant figures. 

This memorandum contains tables giving the values of the right hand sides 
of these equations to six significant figures for the case » = 0. These values have 
been used in an investigation into boundary corrections for axisymmetric cavities 
formed in cylindrical free jets. They are applicable only in the case of a boundary 
which is a free stream function. 

The following tables are also included: 


Table V: —Jo(Ro,my) Table VI: —J1(Ro,my)/Ro, mY 
Table VII: — exp (—o,mx)/Ro, mJ1? (Ro, m) 


where k,, m is the mth positive zero of the Bessel function J,(z). 


1/4[In (r(1 + s)) — 1/2(1 — s)/(1 + s)] — Q(x, », —3)= Cr ne 


A. H. T. 


54[Z].—Proceedings of the First Conference on Training Personnel for the Com- 
puting Machine Field. Edited by Arvid W. Jacobson, Wayne University Press, 
Detroit, 1944. 104 pages. Price $5.00. 
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82 NOTES 


On page 81 of this volume, L. W. Cohen makes two observations. ‘First, the 
effective use of the computing machine depends on the development of appro- 
priate mathematical methods. Second, the development of mathematical methods 
depends on the development of mathematicians.” 

The operators of various computing machines seem to be becoming more 
aware of these dependences, and Wayne University has tried to make information 
available in a series of summer programs. This book reports the second of these 
programs. It is divided into four parts with headings, Manpower Requirements in 
the Computer Field, Educational Programs, Influence of Automatic Computers on 
Technical and General Education, and Cooperative Efforts for Training and Research. 
In the first section an appraisal of manpower requirements in business and 
industry, in government agencies, and by computer manufacturers, is made by 
three contributors. Each of these contributors encourages educational institutions 
to increase their training programs or states his opinion that there will be a tre- 
mendous demand for university graduates trained in the sciences which pertain 
to the computing machine field. 

In the second section of the book several contributors spell out generalities 
of educational programs. In the third part, there is considerably more detail in 
what various universities are doing in their training programs. The Massachusetts 
Institute of Technology lists Machine-aided Analysis, Digital Computer Coding 
and Logic, Numerical Analysis, Methods of Applied Mathematics, a second course 
in Numerical Analysis, Switching Circuits, Digital Computer Applications and 
Practice, Analog Computation, Electronic Computation Laboratory, Pulsed-data 
Systems, and Switching Circuits, as courses directly related to automatic compu- 
tation, for example. 

It would seem well for universities engaged in a training program relating to 
computing machines or contemplating such a program to look through this book 
to extract from the first part some ideas about the potential demand, from the 
second part some ideas about the general aims and content of the program, from 
the third part some ideas concerning appropriate courses, and from the fourth 
part some ideas concerning the possibility of getting help in supporting such 
a program. 

It is unfortunately true that the universities will have a hard time finding 
suitable text material for these courses and capable instructors, for the demands 
reported in the book are consuming the time and the efforts of people who might 
otherwise be available to write texts and to instruct in this field. 


NOTES 


Summary of Educational Opportunities in Electronic Computation 


A tabulated summary of educational opportunities in electronic computation 
is being prepared by Professor H. H. Goode, Professor of Electrical Engineering, 
University of Michigan. This Electronic Computation Education Summary will 
be published in the Transactions of the Professional Group on Electronic Com- 
puters of The Institute of Radio Engineers. 
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NOTES 83 


It is expected that the list will be revised from time to time as conditions 
change. Professor Goode is anxious to make the list as nearly complete as possible, 
and institutions with analog or digital computation courses or facilities, including 
punched card facilities, either planned or in operation for educational purposes, 
should communication directly with Professor Goode. 


Dr. E. W. Cannon Appointed Chief of the Applied Mathematics 
Division of The National Bureau of Standards 


Dr. E. W. Cannon has been appointed Chief of the Applied Mathematics 
Division of the National Bureau of Standards, effective February 28, 1955. 

He will direct the Bureau’s basic mathematical research program which is 
directed principally toward better utilization of electronic computing machinery, 
and progress in numerical analysis, mathematical statistics, and mathematical 
physics. In addition, the Bureau’s Applied Mathematics Division acts as a service 
laboratory to the Bureau and other Government agencies, particularly in the 
fields of numerical computation, statistics, and quality control. 


Policy Committee for Mathematics NBS Technical 
Advisory Committee 


In its report of October 15, 1953, the Ad Hoc Committee, composed of repre- 
sentatives of professional scientific and engineering societies and under the chair- 
manship of M. J. Kelly, appointed by Secretary Weeks to review the activities 
of the National Bureau of Standards, recommended the formation of a set of 
Technical Advisory Committees to advise the Director of the Bureau of Standards 
and his staff on matters which the committees and the staff of the Bureau consider 
of importance. In accordance with this report known as the Kelly report, the 
Policy Committee of the Mathematical Societies of America, which was repre- 
sented on the original Kelly Committee, nominated members for a Technical 
Advisory Committee for the Applied Mathematics Division of the Bureau of 
Standards. 

At present this committee consists of David Blackwell of Howard University, 
Edward U. Condon, consulting physicist, Mark Kac of Cornell University, 
Phillip M. Morse of Massachusetts Institute of Technology, Mina Rees of Hunter 
College (Chairman), and A. H. Taub of the University of Illinois. 

In meetings held on October 23, 1954, and on February 3, 1955, the Advisory 
Committee reviewed the program of the four sections of the Applied Mathematics 
Division, namely the Numerical Analysis Section, the Computation Laboratory, 
the Statistical Engineering Laboratory, and the Mathematical Physics Section. 
It also discussed that part of the work of the Electronic Computer Section having 
a bearing on the work of the Applied Mathematics Division. 

The needs for new computing equipment and methods for financing the acqui- 
sition of such equipment were discussed at both meetings. A survey is being 
made by the Applied Mathematics Division of the computational needs of the 
other divisions of the Bureau of Standards, with two purposes in mind: 1) the 
determination of some of the requirements that the new computing equipment 
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must have, and 2) to acquaint the various divisions of the Bureau of Standards 
with the potentialities of modern computing devices and techniques as re- 
search tools. 

The Applied Mathematics Division in particular and the Bureau of Standards 
in general provide unique services to the scientific and technical community. In 
turn they need the interest and support of this community in facing their financial, 
personnel, and scientific problems. 


The University of Wisconsin Conference 
“The Computing Laboratory in the University” 


The University of Wisconsin is holding a conference entitled “The Computing 
Laboratory in the University,’’ for a two and one half day period beginning 
Wednesday morning, August 17, 1955. There will be a few addresses on the 
computing field in general, several short talks by the ablest users of computing 
equipment, and several panel discussions concerned with the role of computing 
in higher educational institutions. The meeting is planned so as to be of interest 
to administrators and educators in higher educational institutions and to those 
in government agencies and industries who are responsibly concerned with the 
employment of trained personnel. Inquiries concerning the conference may be 
addressed to the Director of the Numerical Analysis Laboratory, 206 North Hall, 
The University of Wisconsin, Madison 6, Wisconsin. 
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CLASSIFICATION OF TABLES 


Arithmetical Tables. Mathematical Constants 
Powers 

Logarithms 

Circular Functions 

Hyperbolic and Exponential Functions 
Theory of Numbers 

Higher Algebra 

Numerical Solution of Equations 


PP 
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Finite Differences. Interpolation 
Summation of Series 

Statistics 

Higher Mathematical Functions 
. Integrals 

Interest and Investment 
Actuarial Science 

Engineering 

Astronomy 

Geodesy 

Physics, Geophysics, Crystallography 
Chemistry 

Navigation 


Aerodynamics, Hydrodynamics, Ballistics 


Calculating Machines and Mechanical Computation 
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